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ABSTRACT 

We present hydrodynamical simulations of the formation, structure and evolution of 
photoionized columns, with parameters based on those observed in the Eagle Nebula. 
On the basis of these simulations we argue that there is no unequivocal evidence that 
the dense neutral clumps at heads of the columns were cores in the pre-existing molec- 
ular cloud. In our simulations, a variety of initial conditions leads to the formation and 
maintenance of near-equilibrium columns. Therefore, it is likely that narrow columns 
will often occur in regions with large-scale inhomogeneities, but that observations of 
such columns can tell us little about the processes by which they formed. The manner 
in which the columns in our simulations develop suggests that their evolution may 
result in extended sequences of radiation-induced star formation. 
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1 INTRODUCTION 

H ii regions are common features of regions of massive star 
formation, such as Orion or the Eagle Nebula, M 16, in which 
ultraviolet photons from young massive stars photoionize 
and photodissociate the molecular gas clouds from which 
the stars originally formed. Emission from the ionized gas, 



in lines! such as the Balmer series or [D IIIJA5UU7A, allows 
the regions to been seen at optical wavelengths, images ot 
these regions have many different forms, although they may 
often be characterised as a photoevaporating blister. 

The basic stratification of gas in an H n region is well un- 
derstood. Close to the central star or stars the gas is highly 
ionized, as a result of the strong photoionizing continuum. 
Atoms which recombine within this region are rapidly pho- 
toionized once again, and the gas is maintained at an equi- 
librium temperature close to 10 4 K. At the edge of this re- 
gion, where the photoionizing continuum spectrum has been 
essentially exhausted, is an ionization front (IF) which sep- 
arates the ionized gas from atomic gas. The ionization front 
is generally narrow, compared to the overall scale of the re- 
gion. Beyond the IF, the gas in the photodissociation region 
(PDR) is kept atomic by photons with wavelengths longer 
than 912 A which do not have sufficient energy to ionize 
hydrogen, but are energetic enough to dissociate molecular 
hydrogen (Bertoldi & Draine 1996). Eventually, the spec- 



trum of photodissociating photons is also absorbed, and the 
gas can remain molecular. 

However, this simple stratification belies the intricate 
structures often observed in H n regions. A feature of par- 



ticular interest is the photoionized columns, or 'Elephant 
Trunks', seen at the edges of regions such as M 16, M 20 and 
NGC 3603. These intrusions of molecular gas into the high- 
pressure environment of the ionized nebula are potentially 
important sites for star formation, and indeed evidence has 
been mounting which associates young stellar sources with 



the tips of the columns in M 16 ( Hester et al. 199q; 



White 



et al. 199S). But the origin, structure and evolution of these 



columns has been the subject of controversy for many years. 
Do they result from pre-existing dense cores in the molecu- 
lar gas or from instabilities in the IF/PDR structure? What 
is the timescale for their evolution? 



1.1 The model of White et al. 



White et al. ( 1999 ) have studied the structure of the columns 
in M 16 using observations in a wide variety of wavebands. 
They suggest that the heads of the columns are pre-existing 
molecular clumps which are now in the process of dynamical 
collapse, driven by an overpressure in the ionized gas around 
the heads, which will soon result in the formation of stars. In 
this model, gas at the heads of the columns in M 16 started 
to collapse at most 10 5 yr ago, when the ionization front (IF) 
moving away from the exciting stars reached their surfaces. 
The velocity of the shock driven into the neutral material 
by the IF is approximately 



pn 



(1) 
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where P sn /kB 



1.2xl0 8 cm- 3 K and P n /k B 



3.5xl0 7 cm _3 K are the pressures inferred for the shocked 
neutral gas and for the gas at the head of the column, re- 
spectively, and p n ~ 2xl0 5 mH 2 cm -3 is the density of the 
gas in the head. White et al. calculate that v B ~ 1.5 km s, 
and so that the shock will cross the head of the column in 
~ 1.5xl0 5 yr. 

In consequence, White et al. argue that the columns 
have been subject to the incident radiation flux for a time 
no longer than this. For instance, the parameters given by 
White et al. imply that the flux of gas into the shock is ~ 
5xl0 10 cm -2 s _1 H nuclei, while the mass flux out through 
the IF is ~ 4xl0 9 cm -2 s _1 H nuclei. Even allowing for the 
different areas of IF and shock, shocked gas should pile up 
between them to a total column of ~ 10 23 cm -2 in 10 5 yr. If 
the clump had been subject to photoevaporation for more 
than 10 5 yr, the shocked gas should now be observable in 
the CO or 850 /im data. 

White et al. provide secondary support for the age they 
infer from detailed chemical and thermal modelling, which 
agrees with the observed properties of the molecular gas, 
and particular the low temperature of the densest molecular 
gas, and from the absence of diagnostics of protostars within 
the cold clumps. They argue that the distance between the 
ionization front and shock propagating into the clump will 
be far smaller than 0.2 pc, so the dense material seen at the 
head of the column cannot itself be shocked gas. 

Based on this model, they propose that the tips of the 
columns in M 16 may represent the best examples known to 
date for earliest stages of Class protostellar development, 
which have not yet collapsed fa r enough to obey the full 
criteria for Class (Andre 1996). 



a) 



b) 



1.2 Other possibilities 



Despite the arguments of White et al. (1999), a lifetime of 
^ 10 5 yr still seems difficult to reconcile with the dynamical 
age of the other structures around the columns. Since the 
density of the ionized gas around the column is far smaller 
than that of the molecular gas, the main ionization front 
is likely to be D-type. As a result, the IF and its leading 
shock will be moving into unperturbed molecular material 
ahead at a velocity ?J lOkms -1 . Indeed, since the age of the 
stellar cluster is 1 — 2 x 10 6 yr, the average velocity of the ion- 
ization front over this lifetime can be no more than 2kms _1 . 
Even if the velocity of the IF were as high as lOkms -1 as 
it propagated along the sides of the column, the IF must 
have passed the head of the column at least 3xl0 4 yr ago, 
so an age substantially smal ler than 10 5 yr seems unlikely. 
Similarly, Hester et al. (1996) argue that the distribution of 



evaporating gaseous globules (EGGs) ahead of the column 
suggests that they are being photoevaporated on a timescale 
of a few x 10 4 yr as the column recedes, so the limit on their 
distribution relies on this photoevaporation timescale rather 
than on the lifetime of the finger itself. Finally, the presence 
of strikingly similar columns in other regions suggests that 
columns may survive for a significant fraction of the total 
age of the H n regions into which they intrude. None of these 
arguments excludes the picture developed by White et al.: 
however, the upper and lower limits on the lifetime are tight, 
and may come to exclude it in time. 

In the present paper, we will explore longer-lived alter- 
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Figure 1. Schematic diagram of the Eagle columns, showing the 
possible positions of the ionization front and shock front. In both 
cases, the ionization front is coincident with the edge of the neu- 
tral material in the column, but in (a), the shock is between the 
ionization front and the dense molecular mate rial in the head of 
the column as disussed by White et al. (1999), while in (b) it is 
between the head and the less dense molecular material in the 
barrel of the column. 



natives to the model of White et al. One is that the dense 
gas in the head is separated from the lower density gas in 
the column by a shock. The comparison between this case 
and the model of White et al. is illustrated schematically 
in Figure hi White et al. rule out such a scenario on the 
basis of the narrowness of the ionization front/shock front 
layer. Indeed, the predicted flux of H nuclei through the 
shock, 10 10 cm -2 s , means that, to have collected the ob- 
served 30 Mq of gas in the head, it would have needed to 
propagate through approximately 1 pc of gas at the mean 
density of the neutral column over a timescale of 4xl0 5 yr: 
substantial fractions of both the size and lifetime of M 16. 
The collapse of the neutral material in the tangential direc- 
tion may, h owever, decrea se the time required to build up 
the mass (cf Bertoldi 198£). Structures of this form may re- 
sult from density structures with initial contrast lower than 
those described by Whi te et al., or fro m partial shadowing 
of the ionization front (Williams 1999). The detailed struc- 
tures found in numerical hydrodynamic simulations below 
are more complex than those shown in Fig. IlL but accord 
broadly with this scenario. 

Another alternative is that the dense gas at the head 
of the column is instead confined by gravity, with the 
shock leading the ionization front having effectively by- 
passed the region and moved away downstream. White 
et al. find that self-gravity may be important at the 
head of the column. The size of the region dominated 
by the gravity of the observed molecular core is about 



0.13(M/30M Q )(cr/lkms~ 1 )" 2 pc, where a ~ lkms -1 is 
the effective sound speed of the gas (including contributions 
from turbulence and magnetic fields as well as the thermal 
sound speed) and M ~ 30 Mq is the mass of neutral gas in 
the head of the column. 

In the following sections, we discuss numerical hydro- 
dynamic simulations of these scenarios. These have been tai- 
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lored to the central column in M 16, which has the most sym- 
metrical appearance and for which the observational data 
have the highest spatial resolution. We follow Hester et al.'s 
nomenclature, and term it column II (it is referred to as II2 
by White et al.). 

The structures which we find have some similarities 
to the equilibrium come tary cloud solutions discussed by 
Bertoldi & McKee (199C), based on the collapse of isolated, 



near-spherical clumps. However, our time-dependent hydro- 
dynamical models give useful information on flow stability 
and the presence of shocks. Simil ar str uctures have been 
found by Garcfa-Segura & Franco (1996), in their study of 



the development of H 11 regions in inhomogeneous environ- 
ments, although as a result of the global scope of their sim- 
ulatio ns ind ividual columns were poorly resolved. Mellema 
et al. ( 199£ ) model the formation of neutral tails behind iso- 
lated dense clumps in some detail. The models we present 
here have been developed for cases where the columns are 
close to the edge of an H 11 region, with parameters chosen 
to compare directly with those in M 16. 

In Section H, we discuss the derivation of suitable pa- 
rameters for our numerical models from the observational 
data, and particularly the comparison of the pressures of 
ionized and neutral gas across the ionization front at the 
head of the column. In particular, we find that the observa- 
tional evidence for an overpressure around the heads is weak. 
In Section H, we discuss the numerical method we have used, 
and, in Section W, the initial conditions for the simulations. 
In Section H, we present the results of these simulations, and 
in Section discuss their implications. 



2 THE COLUMNS IN THE EAGLE NEBULA, 
M 16 

2.1 Observational data 

The columns in M 16 have been studied by many authors. 
Perhaps the best op tical data is the HST narrow band imag- 
ing of Hester et al. (1996). We retrieved two Ha images from 
this data set from the HST archive, and stacked and median- 
limited them to remove strong cosmic ray features and very 
fine-scale structures. The resulting image is shown in Fig- 
ure H. There are many striking irregularities in this image. 
Small clumps break the surface, and some appear fully sepa- 
rated from it. Close to the apex, the surface is puckered on a 
fine scale. Some of the clumps harbour K-band sources, and 
are identified as protostars by Hester et al. Away from the 
head, the surface irregularities seem to grow in width, with 
bright regions where their surfaces face towards the source of 
ionization. In addition, we note the low-magnitude striations 
which are seen along lines approximately perpendicular to 
the nearest surface. 



White et al. (1999) combined this optical data with ob- 



servations of M 16 in molecular lines, millimetre and sub- 
millimetre continuum and the mid-infrared, to present a re- 
markably complete picture of properties of the gas and dust 
in and around the columns. Within the columns, molecu- 
lar gas is probed by line emission and the associated dust 
is seen in millimetre/sub-mm emission, while Ha and radio 
continuum observations are sensitive to structures in the 
surrounding ionized gas. 




Figure 2. Cleaned Ha image, shown negative and with a 
grayscale saturated in order to bring out low contrast features. 
The axes are labelled in pixels, where each pixel corresponds to 
1.4xl0 15 cm, for Hester et al.'s adopted distance of 2kpc. 



Levenson et al. (2000) have observed excited H2 emis- 
sion in M 16, to study the properties of the PDR which 
leads the ionization front. They find that their observations 
are consistent with a stationary photodissociation region, 
around 2xlO lr cm deep (although since the properties of 
the incident radiation field are not well determined, it is 
possible that shocks driven by the IF might also contribute 
to the collisional component of this emission). 



2.2 Comparison with photoevaporated winds 

In Figure H, we show cuts through this image. The first is 
along the axis of the large-scale column which is vertical 
in Fig. bl The values are averages for a strip between 345 
and 355 on the i-axis. The second is perpendicular to this, 
averaged for a strip between 395 and 405 on the y-axis. 

In Fig. kl we show the radial emission measure profile 
derived from i sother mal wind models, as described by Hen- 
ney & Arthur (1998). These authors found that such models 
gave a very good fit to Ha emission profiles of gas photoevap- 
orated from proplyds in the Orion nebula. We have here been 
rather less strenuous in the correction of the observational 
data for contamination by continuum and [N 11] emission, 
but our aim in this section is simply to get an indication of 
the properties which pertain at the head of column II. 

Comparing the model profiles with Fig. H, we see that 
their overall forms are very similar. The central brightness 
of the columns appears rather lower, relative to that at large 
distances, than the model would predict, but that may be 
explained by the decreasing intensity of illumination as the 
surface of the column points away from the head and by 
the occulting of background emission by the column. The 
sharpness of the inner edge of the profile suggests that the 
column is slightly closer to us than is the photoionizing star. 

More fundamentally, if (following Hester et al.) we take 
the radius of curvature of the head of the column to be 
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Figure 3. Sections through the data (a) top-to-bottom (b) left- 
to-right. 




Figure 4. Sections through spherical flow models for isothermal 
winds which start from Mach 1, 1.2 and 2 at the surface of the 
column, and for constant velocity flow. The plots have been nor- 
malized to fixed peak emission measure: the sharpest plot is for 
a Mach 1 outflow while the smoothest is for constant velocity 
outflow. 



90 pixels, then the breadth of the observed Ha profile is 
difficult to account for. Not even the most extreme of the 
models in Figure W explains the width of the tran sverse 
light distribution. For comparison, Sankrit & Hester (200C) 
find that power-law models with a radius of curvature of 
2xl0 17 cm ~ 150 pixels fit the emission profile in sev- 
eral emission lines, which confirms that the effective radius 
of curvature is larg er th an that apparent in the images. 
Megeath & Wilson (1997) also find that the radio emission 



measure is more extended in their observations of photoe- 
vaporating clumps in NGC 281 West than for a n e oc 1/r 2 
density model. 

One possible explanation for the M 16 results is that the 
head of the observed column is more elongated along the line 
of sight than perpendicular to it. Alt ernatively, if the column 
is directed into the plane of the sky (Pound 199£), the inner 
peak of the emission measure profile may be hidden from 
view, although the inclination does not appear to be that 
extreme. 



2.3 Inferred density and radiation field 

For the model ionized flow profiles, isothermal flow with ini- 
tial Mach numbers of 1, 1.2, 2 and constant velocity flow, the 
effective integration distance (i.e. the peak emission measure 
divided by the square of the peak density) is 0.8, 1.0, 1.3, 
and 7r/2 times the radius of curvature, respectively. This 
may be compared to the value of one-half assumed by Hes- 
ter et al. (since the profiles are well resolved, the correction 
due to resolution effects is minimal) . When the uncertainty 
in the local radius of curvature is included, it seems very 
likely that Hester et al.'s value of 4000 cm -3 for the electron 
density at the head of column II may have been overes- 
timated. Indeed, Hester et al.'s own photoionization mod- 
elling suggested that the density at the head of the column 
is 2000 cm' 3 (although th is has subsequently been revised 
by jankrit fc Hester 2000 ). As a result, the ionized gas may 
not in fact have a greater pressure than the effective pressure 
in the neutral gas at the head of the column. 

The density at the head of the column can also be 
used to determine the incident ionizing flux. The thickness 
parallel to the radiation field is roughly I ~ 0.04 pc from 
Figure H. Balancing the number of recombinations in the 
flow with the intensity of the incident ionizing radiation, J, 



suggests that J ~ aBn„I ~ 1.3X10 11 cm _2 s _1 , for a peak 
electron density of n c ~ 2000 cm -3 . It is difficult to fol- 
low White et al.'s analysis of their VLA data, but assuming 
that the integrated 48 mjy emission corresponds to that oc- 
culted by the head of the column, wh ich has a radius of 
curvature of 0.04 pc (Hester et al. 1996), the incident ioniz- 



ing intensit y pred icted by the formulae of Lefloch, Lazareff 
& Castets ( |1997D is J ~ 1.5xlO n cm' 2 s~\ These values 
are mutually consistent, but smaller than those found by 
White et al., and also smaller than the values, inferred from 



stellar censuses, of 4.4x10 cm s (Hester et al. 199(:) 



or 2.6x10 cm s (from the total flux given by 



Pound 



199£). However, this difference can easily be accounted for 
by uncertainties in the ionizing flux of the stars, as the result 
of dust absorption or photoionization between the ionizing 
stars and the columns, or by the geometry of the region. 

In our models, we neglect the role which a PDR may 
have in the dynamics of the region. The depth of the PDR 
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in M 16, as derived by Levenson et al. (2000), is close to 



the radius of curvature of the finger tips, so this assumption 
is only marginally valid. However, since the column density 
inferred is close to the stationary models (rather than being 



cf. 



Bertoldi & Draine 



subs tan tially affected by advection 

199 6J), [ t seems likely that the jump conditions across the 
combined IF/PDR structure are determined principally by 
the flux incident on the IF.r] 



2.4 Derived parameters 

Based on the observations, we will adopt the following 
parameters for our modelling of column II. The pressure 
in the ionized gas at the head of the clump is Pi/fce — 
3xl0 7 cm~ 3 K; this ionized gas is assumed to have an 
isothermal sound speed of lOkms" 1 . The density in the 
molecular gas in the core at the head is n(H.2) = 2x 10 5 cm~ 3 
while that in the column of gas behind is typically 10 times 
smaller than this. White et al. suggest that while the tem- 
perature of the gas is 20 K throughout the column, the ef- 
fective pressure within the columns is dominated by the ef- 
fects of small-scale magnetic fields and turbulence, as evi- 
denced by the smooth profiles of molecular lines with widths 
~ 2kms _1 . We cater for this in our dynamical models by 
assuming that the neutral gas behaves as if it were isother- 
mal with a temperature of 200 K (note that the turbulent 
contribution to the pressure in the molecular gas also re- 
duces the importance of heating in the PDR for dynamical 
modelling). The total mass in the core of column II is 31 Mq 
and its radius is 0.085 pc. We assume a photoionizing flux of 
10 11 cm -2 s _1 , impinging on the column parallel to its axis. 



3 NUMERICAL METHOD AND 
ASSUMPTIONS 

We u se the hydrodynamical code described by Williams 
( 1999 ), which is based on the second-order Godunov al- 



gorithm by Falle (1991). Ionization and recombination are 



treated by including the equations 



da; 

''dt 

dJ 

dz 



an(l — x)J — aan x 
— an(l — x)J. 



(2) 
(3) 



Here n is the density of hydrogen nucleons, x the ioniza- 
tion fraction and J is the number flux of ionizing photons, 
assumed to propagate parallel to the z axis. The photoion- 
ization cross-section of hydrogen, a, and the case B recombi- 
nation coefficient, ob, are taken as constants. The electron 
density, nx, is included as an advected variable in the hydro- 
dynamic solver, and the ionization and recombination equa- 
tions included by operator-splitting, implemented at the end 
of each step (and also for the intervening half-steps). The 

* If, in other cases, the depth of a PDR in pressure equilibrium 
with the photoionized flow were larger than the radius of a glob- 
ule, the PDR would be important in development of the globule; if 
it were larger than the inter-globule separation, the development 
of photoionized columns might be suppressed, with the internal 
structure of the molecular cloud only observable in atomic and 
molecular species. 



equations are integrated using an implicit scheme to advance 
x in time and J across each grid cell. 

The code does not include a treatment of the diffuse 
ionizing radiation field (beyond using the case-B recombi- 
nation coefficient). In the on-the-spot approximation, the 
diffuse field is found to be ~ 15 per cent of the local direct 
radiation field (Canto et al. 199SJ). The characteristic length 



for reabsorption of diffuse photons is I ~ J/asn ~ 0.05 pc 
(dependent, or course, on spatial variations of J and n), so 
the diffuse field might be expected to have some effect even 
on structures as large as the columns. More detailed treat- 
ments predict that the diffuse radiation is a smaller fraction 
of the direct field away fr om the edge of the global photoion- 



ized region ( Rubin 1968 ). Close to the edge, where relative 
size of the diffuse field is larger, the soft spectrum of the dif- 
fuse field will result in a smaller path length to absorption 
than for the direct field. The low contrast of the surfaces 
of the observed columns which do not point towards the 
ionizing source is empirical evidence that the diffuse radi- 
ation field would is adequately treated by the on-the-spot 
approximation in the present case. 

Since the internal structures of IFs are not resolved in 
the simulations, most of the gas will be either almost com- 
pletely ionized or fully neutral. Intervening regions with in- 
termediate ionization are likely to result principally from nu- 
merical mixing, so physically detailed schemes which assume 
the gas within the computational cells is reasonably uniform 
are not appropriate. Instead, we set isothermal sound speed 
in the gas to 



= (l + 99z) 1/2 kms~ 



(4) 



This expression results from assuming that in cells with 
< x < 1, the gas is split into separate regions of fully 
ionized and fully neutral gas, in pressure equilibrium with 
each other, and so gives a better account of the sub-grid 
structure at the (dynamically significant) D-type IF. 



4 INITIAL AND BOUNDARY CONDITIONS 

The main model presented here is calculated in cylindrical 
symmetry, with coordinates (r, z). A uniform 128 x 448 grid 
covers physical dimensions < r < 0.5 pc and < z < 
1.75 pc. The boundary conditions are reflective on the axis, 
r = 0, as required by symmetry. On the boundary towards 
the source, z = 1.75 pc, free-flow boundary conditions are 
used when the flow is outwards, mirror conditions when it 
is inwards, although as the outward flow is generally super- 
sonic this has little influence on the flow inside the grid. 
Along the outer radial boundary, r = 0.5 pc, the conditions 
beyond the grid are either free-outflow/zero-inflow if the ion- 
ization fraction is larger than 0.1, or mirror conditions if 
it is less than this value. This somewhat artificial criterion 
changes the behaviour at the boundary at the IF, in order 
to model an ionized flow which is divergent on large scales. 
The effect depends little on the specific limiting value, as 
away from the IF the gas is in general either almost entirely 
ionized or almost entirely neutral. 

On the edge of the grid towards the molecular cloud, 
z — 0, the boundary condition is set by assuming that the 
material retains its initial state, i.e. in most of our exam- 
ples a flow of gas towards the IF. As the solution develops, 
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a shock often becomes captured at the boundary, which re- 
duces the inflow speed and increases the inflow density. The 
flux of neutral gas into the grid through the plane z = is, 
in all cases, small enough that it could not prevent the IF 
from entering if the entire grid were full of neutral material. 
In many cases, the density of the ionized gas becomes large 
enough to absorb all the ionizing photons within the size of 
the grid. Where the 'molecular' boundary covers the entire 
upstream edge of the grid, a plane-parallel IF is an equilib- 
rium solution. If the ionized flow is allowed to escape from 
the side of the region, it is not obvious what the form of the 
equilibrium will be, if indeed there is one. One possibility is 
that the IF surface will remain smooth, and a steady solu- 
tion will be found in which the variation in pressure at the 
IF is balanced by changes in flow velocity via the Bernoulli 
effect. 

However, these IF are believed to be unstable to cor- 
rugation modes on both small and large scales, although 
for wavelengths larger than 6xl0 15 nJ 1 cm recombination 
decrease s the grow t h rate of these instabiliti es decreases 
rapidly ([Kahn 1958J; |Axford 1964} |Sysoev 1997|). This wave- 



length is comparable to our grid scale of 10 cm for typi- 
cal ionized gas densities at the IF of ~ 10 3 cm -3 . We now 
present numerical models which illustrate what may happen 
for a range of initial conditions. 

The models we have calculated have initial conditions 
split into two regions: ionized gas with a density n(H total) = 
200 cm~ 3 towards the source of ionizing photons and moving 
towards it at lOkms" 1 relative to the neutral gas behind. 
This neutral gas has density, n(H total) = 2x 10 4 cm -3 (n.b., 
not n(H2), see below). Except where noted, in the initial 
conditions the neutral gas is assumed to move at 2kms _1 
towards the ionizing source in order to keep the ionization 
front broadly at rest within the grid, with the ionized gas 
assumed to stream towards the radiation source lOkms -1 
more rapidly than this. Pressure equilibrium holds where 
ionized and neutral gas are alongside each other. We take the 
ionizing flux incident at the top of the grid in the figures to 
be 1.5x 10 11 cm -2 s _1 , since recombinations in the upstream 
flow are found to reduce the ionizing radiation incident on 
the head of the column somewhat. The Stromgren distance 
in the ionized gas, ~ 3pc, is sufficiently great that a column 
the length of the computational grid can be maintained with 
almost complete ionization for the initial conditions. 

The boundary between these regions is assumed either 
to be plane (to study the development for small perturba- 
tions) or to include a column of dense gas (to study the 
likely long-timescale development of the columns observed 
in many regions). 



5 RESULTS 

We will now present the results of four different scenarios 
for the development of photoevaporated columns. 



Case I is based on the model of White et al. (1999), in- 
cluding initial density stratification and gravity, 

Case II has an initial short cylindrical column, 

Case III has a 10 per cent deficit in photoionizing flux close 
to the axis, 

Case IV has dense neutral gas initially filling the grid, with 
gravitational forces from a 30 Mq object included. 



Case I: In Figure H, we present the results of a first 
simulation in which the flow conditions modelled on those 
inferred by White et al. The initial conditions have a dense 
core of material at the head, and a lower density col- 
umn behind. A gravitational potential from a 30 Mq mass, 
smoothed at a radius of 0.08 pc, is centred on the core of the 
head. The uniform gas at the head of the column initially 
falls in, as it is not in equilibrium with the gravitational 
field, but does not collapse far because of the smoothed core 
of the field. 

By the time of the second frame, 10 5 yr after the start of 
the simulation, the shock at the head of the the column has 
passed through much of the dense gas. The shock proceeds 
down the column in later frames. Soon the ram pressure of 
the radiation field is sufficient to push the dense gas away 
from the gravitating core. The amount of gas with density 
above 10 5 cm' 3 in these models remains roughly constant, 
and so the gravitational field would in reality respond to this 
movement. Nevertheless, it is striking that the density field 
in the column remains close to that inferred by White et al. 
even when the head is a considerable distance away from 
the gravitating centre, rather than requiring the head of the 
column to be caught within 10 5 yr after the IF reaches it. 

In the remaining examples, we investigate possibilities 
by which such structures may have formed, and what con- 
ditions are necessary for them to be long-lived. 



Case II: In Figure pi we illustrate a case in which a 
short column of gas is initially present. The radius of the 
column is initially 0.125 pc and it stands proud by 0.5 pc 
ahead of the main body of molecular material. The column 
is again seen to be long lived, but highly variable. Certain 
broad features characterise the structure. The highest den- 
sities in the neutral gas are seen close to the head of the 
column, and reach values similar to those inferred for the 
column tip by White et al. 

The surface of the column is convoluted. Neutral clumps 
split off from from the column at different stages (e.g. close 
to the base of the column in Figure of) and linear struc- 
tures are seen where different streams of ionized gas interact. 
These are reminiscent of the observed clumps and striations 
in the flow (see Figure H and Section H). The striations are 
sometimes seen to be triggered by surface inhomogeneities 
resulting from the discrete nature of the grid. However, the 
growth of these instabilities appears to be a generic feature 
of the development of obliquely-illuminated ionization fronts 
(Vandervoort 1962), even when account is taken of the ef- 
fects of recombination in the upstream flow (Williams, in 
preparation). 

In the frames shown in Figure pi the dense column be- 
comes more marked with time. Initially, the ionizing field 
generates a strong flow away from the head of the column, 
a low density region ahead of it and a net acceleration away 
from the ionizing source (the well-known 'rocket effect', Oort 



Spitzer 1954). However, we see here that short lengthscale 



which characterises the photoevaporative flow from the head 
means that its pressure decreases rapidly, and eventually a 
termination shock forms. The convergence of this termina- 
tion shock at the axis leads to the formation of an over-dense 
region at larger distances from the column. Such over-dense 
regions lead to the formation of some thin filaments con- 
nected to the main ionization front [e.g., those seen close to 
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(b) 
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f/pc 



r/pc 



Figure 5. Greyscale plots of log density for Case I: an initial column of gas based on the model of White et al. The plots are shown at 
the labelled times after the start of the simulation. The density range calibration is shown in the colour bar; a single contour marks the 
location of the ionization front. 



the base of the column in Figure |g(b)] . Filaments at the head 
of the column can allow the column to move forward once 
again, although the rigidity provided for the filaments by 
the reflecting boundary condition on-axis may have a signif- 
icant influence by preventing non-axisymmetric instability 
modes. 

Case III: In Figure M, we show the development of 
columns for a case where there is a 10 per cent weakening 



of the impinging ionizing radiation field for a small radius 
around the axis. A narrow feature in the radiation field such 
as this might formed by an upstream overdensity: a photo- 
evaporating globule, for instance, or dense structures close 
to the stars caused by interacting stellar winds. 

A strong column soon forms in this case. Once it has 
formed, its evolution is not greatly influenced by the small 
perturbation in the radiation field. The structure is broadly 
similar to those seen in the previous two cases, although the 
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Figure 6. Grcyscale plots of log density for Case II. The plots are shown at (a) initial time (b) 10 5 yr, (c) 2xl0 5 yr, (d) 3xl0 5 yr, (e) 
4xl0 5 yr, and (f) 5xl0 5 yr after the start of the simulation. The full density range is typically from 3 to 10 6 cm - 3 , although the gas 
which fills most of the volume of the column is between 5x 10 4 and 5x 10 5 cm -3 . 



column has a somewhat higher aspect ratio and the shock 
remains within the computational grid for a longer period. 



Case IV: In Figure we take uniform dense gas to 
nearly fill the grid, and include a gravitational field equiv- 
alent to a 30- Mq core smoothed at a radius of 0.02 pc. In 
this case only, the neutral gas is initially at rest. As the 
simulation develops, mass first falls into the gravitational 



potential, to form a dense core close to hydrostatic equilib- 
rium. A shock and ionization front driven by the imping- 
ing ionizing radiation propagate over this core, generating 
a column structure. This column is very long-lived, tend- 
ing to form a cometary structure at late times. The length 
of this cometary structure may in reality be limited by the 
non-plane parallel ionizing radiation field resulting from dif- 
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(a) 



(b) 



(c) 




Figure 7. Greyscale plots of log density for Case III, stationary, initially uniform gas with 10 per cent decrease in ionizing radiation field 
within a 0.125 pc radius of the axis. The plots are shown at (a) initial conditions, (b) 1x10 s yr, (c) 2 x10 s yr, (d) 3 x10 s yr, (e) 4 x10 s yr, 
and (f) 5 x10 s yr after the start of the simulation. The density range calibration is shown in the colour bar. 



fuse em ission and the multiple ionizing stars ( Pavlakis et al 
200l|)~l 



ways present within the solution, but are a characteristic of 
the dynamically disturbed nature of the flow. 



This model seems similar to that discussed by 
White et al. However, it will be noted that by the time the 
ionization front reaches the dense clump, it interacts with 
the reflection of the main shock front from the dense core, 
rather than driving a separate shock itself. Shocks are al- 



Here we have imposed a separate potential, rather than 
calculating the gravity self-consistently, primarily as a com- 
putational convenience. However, it should be noted that 
in the Eagle columns, a reasonable fraction of the gravita- 
tional potential may be generated by protostellar clumps 
which have effectively decoupled from the flow. The small 
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(a) 



(b) 



(c) 




Figure 8. Greyscale plots of log density for Case IV: stationary, initially uniform gas with a 30 Mq gravitational field applied. The plots 
are shown at (a) 2 x10 s yr, (b) 3xl0 5 yr, (c) 4xl0 5 yr, (d) 5xl0 5 yr, (e) 6xl0 5 yr, and (f) 7xl0 5 yr after the start of the simulation. The 
density range calibration is shown in the colour bar. 



smoothing radius means that the gravitational field could 
not be generated by the self-gravity of a stable isothermal 
molecular gas core. It is possible that protostellar cores may 
be formed from gravitational instabilities in the swept up gas 
shell (Francis & Whitworth, in preparation), or be present 
in the initial molecular cloud. 



gravity, then these cores must have been surrounded by sig- 
nificantly hotter, less dense gas if the overall structure was 
not to be unstable to gravitational collapse. In this case, 
the source of the ~ 10 4 cm -3 gas in the column behind the 
column is unclear. 



Note, however, that if the heads of the columns were 
formed from individual isothermal cores confined by self- 
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r/pc 

Figure 9. Greyscale plot of Ha emission. The values are calcu- 
lated for the simulation frame shown in Figure 0(f), 5 X 10 5 yr after 
the start of the simulation. The greyscale is linear, and saturates 
at an intermediate flux. Note that the emissivity in the symmetry 
plane is shown here, rather than the line-of-sight integral. 



shows the presence of bright rims and strong linear features 
in the flow away from the surface. 

Since the structures appear to be close to equilib- 
rium, the observations contain limited information about 
the mechanisms by which they formed: as we have shown, 
a variety of initial conditions can lead to very similar final 
structures. However, the observations seem not to require 
highly stratified density structures to have been present 
in the original molecular cloud. Indeed, it seems possible 
that the columns result from instabilities of the rim of a 



champagne-type expandi ng H n region ( Tenorio-Tagle 197!: ; 
Bodenheimer et al. 1979 ). It is difficult to make substantial 



analytic progress in the study of such surface perturbations, 
as a result of the nonlocal nature of the system. As we have 
seen in the simulations, density enhancements resulting from 
fluctuations in the surface a significant time ago can result in 
long term variations in the ionizing flux incident on the sur- 
face a considerable distance from the original perturbation. 
While it is certainly true that the assumption of cylindri- 
cal symmetry enhances the variation of ionizing flux on the 
axis, it seems likely that the long lags between perturbation 
and response inherent in the system will in general lead to 
overstabilities of the flows, when studied on large scales. 

The interpretation we have developed, however, has its 
own coincidental features. Most particularly, the mass of 
shocked molecular gas at the head of the columns is curi- 
ously close to the Jeans mass. There is no obvious dynamical 
reason for the mass of dense gas to settle on any particu- 
lar value: indeed, we have seen the mass of dense gas to 
be highly variable in our simulations. A possible explana- 
tion for why the mass is close to the limit of collapse is 
that, in fact, the mass at the head increases as the column 
erodes, but once it reaches the limit of stability it collapses. 
The dense, self-gravitating globules produced will then no 
longer be subject to significant radiation or hydrodynamic 
forces, and can be left behind by the further photoevapo- 
ration of the column. In this manner, an ionization front 
can lead to stars forming in an episodic manner over a sub- 
stantial period, compared to the one-off radiation-induced 
collapse of clumps which has been studied previous l y (Bok, 



5.1 Discussion 



as referred to by Port 1954 ; gandford et al. 1982 ; Bedijn 
Tenorio-Tagle 1984| ; |Bertoldi 1989[ ). Note that the ion 



In each of the cases presented in the present section, struc- 
tures strikin gly reminiscent of those observed by White 
et al. (1999) are maintained within the computational grid 
for periods of several 10 5 yr. The overall structures observed 
for models run for a range of similar parameters were hroadly 



izing radiation incident on a clump will increase gradually 
when the main ionization front is of D-type, rather than 
suddenly as often assumed. The shock which precedes the 
IF may however r esult in a thin gas layer unstable to gravi- 



o(h 
anj. 



the sa me, suggestion that our results arc not dependent on 



tational collapse (Elmegreen & Lada 1977; Whitworth et al 



the details of our numerical or physical assumptions. This 
suggests that structures seen in M 16 are in a near-steady 
state, with the head of the column confined by the rocket 
effect at the front and by the ram pressure of the gas be- 
ing swept up, by self-gravity, or a combination of the two. 
The structures can survive for several times 10 5 yr, longer 
than the shock crossing time of the head of the column, and 
similar to the dynamical age of M 16 as a whole. However, 
the influence of the symmetry axis on the stability of these 
columns means that this interpretation is tentative, until 
fully three-dimensional models are available. 

In Figure |9_[ we show the Ha emissivity expected from 
the gas within Case I after 5 x 10 5 yr fro m the start of the 
simulation, calculated as in Steffen et al. (1997). This image 



1994) and may trigg er the collapse of any pre -existing self- 
gravitating clumps (Megeath & Wilson 1997). To develop 



these ideas further requires models in which the self-gravity 
of the molecular gas is included. 



6 CONCLUSIONS 

In this paper, we have presented hydrodynamical models of 
the development of photoionized columns, which agree well 
with the observations of the columns in M 16. Our mod- 
els contrast with the results of White et al., who interpret 
observations of the columns in M 16 as requiring that they 
have been exposed to the radiation field for less than 10 5 
years, the dense material in the head of the column being a 
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gravitationally-bound condensation from the initial molecu- 
lar cloud. 

It remains to determine unambiguously how columns 
such as those in M 16 form, and how they evolve, if indeed 
there is any single mechanism. It seems likely, however, that 
these structures are often long lived components of the en- 
vironment of massive stars, and of great importance in the 
formation of lower mass stars in these environments, and so 
they warrant further detailed study. 



Williams R.J.R., 1999, MNRAS, 310, 789 
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